source('00_util_scripts/mod_bulk.R')
source('00_util_scripts/mod_bplot.R')
library(readxl)

excel_sheets('mission/SLE_TRPM2_MfMo/data/GSE255108_edgeR_DEGs_2021-03-24.xlsx')

lps_vs_dmso_mm <-
read_excel('mission/SLE_TRPM2_MfMo/data/GSE255108_edgeR_DEGs_2021-03-24.xlsx', 4)

lps_vs_dmso_mm |>
  filter(gene_name %in% gene_of_int) |>
  select(gene_name, contains('LPS')) |>
  pivot_longer(-1) |>
  mutate(type = str_extract(name, 'logFC|PValue|FDR'),
         group = str_extract(name, '\\d+h')) |>
  filter(type != 'PValue', !str_detect(name, 'YKL')) |>
  pivot_wider(id_cols = c(gene_name, group),
              values_from = value, names_from = type)

gene_of_int <- c('Trpm2','Ccr2','Hmgb2','Cd14','Tnf')

lps_vs_dmso_cpm <-
  read_excel('mission/SLE_TRPM2_MfMo/data/GSE255108_edgeR_DEGs_2021-03-24.xlsx', 3)

lps_vs_dmso_cpm |>
  filter(gene_name %in% gene_of_int) |>
  select(gene_name, matches('^LPS|^DMSO')) |>
  pivot_longer(-1) |>
  mutate(group = str_extract(name, '\\d+h') |> fct_relevel('12h',after = Inf),
         value = log1p(2^value)) |>
  ggplot(aes(group, value, fill = group)) +
  stat_mean(geom = 'col') +
  geom_jitter(width = .1, height = 0) +
  facet_wrap(~gene_name, scales = 'free_y') +
  labs(title = 'Mouse BMDM after LPS treatment', y = 'log(tpm+1)',
       subtitle = 'GSE255108') +
  stat_compare_means(ref.group = '0h', label = 'p.signif',
                     method = 't.test') +
  theme_pubr() +
  scale_y_continuous(expand = expansion(mult = c(NULL,.2)))

# GSE286554 4h LPS + 30m ATP ------------
lpsatp <-
read_delim('mission/SLE_TRPM2_MfMo/data/GSE286554_gene_sample_FPKM_with_symbol.txt.gz')

lpsatp |>
  filter(Symbol %in% gene_of_int) |>
  select(-c(1,3)) |>
  pivot_longer(-1) |>
  mutate(group = ifelse(str_detect(name, 'AT'), 'LPS+ATP', 'PBS') |>
           fct_relevel('PBS'),
         genotype = ifelse(str_detect(name, 'WT'), 'WT', 'Nlrp3-KO') |>
           fct_relevel('WT'),
         Symbol = fct_relevel(Symbol, 'Trpm2','Ccr2','Hmgb2','Tnf')) |>
  filter(genotype == 'WT') |>
  ggplot(aes(group, log1p(value), fill = group)) +
  stat_mean(geom = 'col', position = 'dodge2') +
  geom_jitter(height = 0, width = .1) +
  facet_wrap(~Symbol, scales = 'free_y') +
  labs(title = 'WT BMDM treated with 4h LPS followed 30min ATP',
       y = 'log(fpkm+1)', subtitle = 'GSE286554') +
  theme_pubr() +
  scale_fill_hue(direction = -1)

# GSE272546 ATP 3h -------------
GEOquery::getGEOSuppFiles('GSE272546', makeDirectory = F, fetch_files = F) |>
  pull(url)

thp1_atp <- readxl::read_excel('GSE272546_THP-1_macrophages_processedData_ATP...Control.xlsx')

thp1_atp |>
  filter(SYMBOL %in% str_to_upper(gene_of_int)) |>
  ggplot(aes(SYMBOL, log2FoldChange, fill = padj)) +
  geom_col()

dummy_count <- thp1_atp |>
  filter(SYMBOL %in% str_to_upper(gene_of_int)) |>
  select(baseMean, log2FoldChange, lfcSE, padj, SYMBOL) |>
  mutate(fold_to_pbs = 2^log2FoldChange,
         pbs_mean = baseMean * 2 / (1+fold_to_pbs),
         atp_mean = pbs_mean * fold_to_pbs) |>
  rowwise() |>
  mutate(pbs_value = list(rnorm(n = 3, mean = pbs_mean, sd = lfcSE * 6)),
         atp_value = list(rnorm(n = 3, mean = atp_mean, sd = lfcSE * 6))) |>
  select(SYMBOL, pbs_value, atp_value) |>
  unnest(c(pbs_value, atp_value))

dummy_count |>
  pivot_longer(-1) |>
  mutate(group = str_extract(name, 'pbs|atp') |> str_to_upper() |>
           fct_relevel('PBS')) |>
  ggplot(aes(group, log1p(value), fill = group)) +
  stat_mean(geom = 'col') +
  geom_jitter(width = .2, height = 0) +
  facet_wrap(~SYMBOL, ncol = 4) +
  theme_pubr() +
  labs(title = 'THP-1 cells treated with 3h ATP', y = 'log(tpm+1)',
       subtitle = 'GSE272546') +
  scale_fill_hue(direction = -1)

# GSE260996 BMDM cocktail ----------
GEOquery::getGEOSuppFiles('GSE260996', F, fetch_files = F)

bar <- 'mission/SLE_TRPM2_MfMo/data/bmdm_cocktail/GSM8129359_CF1Aligned.RData' |>
  read_rda()

ent2sym <- cf1$annotation$GeneID |>
  clusterProfiler::bitr(fromType = 'ENTREZID', toType = 'SYMBOL',
                        OrgDb = 'org.Mm.eg.db')

ent2sym |> head()

cf1$counts |>
  as_tibble(rownames = 'ENTREZID')

cctl_path <-
  list.files('mission/SLE_TRPM2_MfMo/data/bmdm_cocktail', full.names = T)

cctl <- cctl_path |>
  str_extract('.{3}(?=Align)') |>
  set_names(cctl_path, nm = _) |>
  map(read_rda, .progress = T)

cctl_df <- cctl |>
  map(\(x)x$counts |> as.data.frame()) |>
  list_cbind()

cctl_long <- cctl_df |>
  as_tibble(rownames = 'ENTREZID') |>
  unnest(-1) |>
  right_join(ent2sym) |>
  relocate(SYMBOL) |>
  select(-ENTREZID) |>
  pivot_longer(-1)

cctl_long <- cctl_long |>
  mutate(name = str_extract(name, '^...'),
         group = str_extract(name, '^..') |>
           case_match('CF' ~ 'CpG+PBS',
                      'CL' ~ 'CpG+LPS',
                      'CT' ~ 'CpG+TNF',
                      'FC' ~ 'PBS+CpG',
                      'FF' ~ 'PBS+PBS',
                      'FL' ~ 'PBS+LPS',
                      'FP' ~ 'PBS+polyIC',
                      'FT' ~ 'PBS+TNF',
                      'PF' ~ 'polyIC+PBS',
                      'PL' ~ 'polyIC+LPS',
                      'PT' ~ 'polyIC+TNF'))

cctl_long |>
  calc_tpm(sample = name, abundance = value) |>
  filter(SYMBOL %in% gene_of_int,
         str_detect(group, '^PBS')) |>
  mutate(group = str_remove(group, 'PBS\\+') |> fct_relevel('PBS'),
         SYMBOL = fct_relevel(SYMBOL, 'Trpm2','Ccr2','Hmgb2','Tnf')) |> 
  ggplot(aes(group, log1p(tpm), fill = group)) +
  stat_mean(geom = 'col') +
  geom_jitter(height = 0, width = .2) +
  facet_wrap(~SYMBOL, scales = 'free_y') +
  theme_pubr(x.text.angle = 90) +
  labs(title = 'm-CSF conditioned BMDM treated with 4h ligands',
       subtitle = 'GSE260996', y = 'log(tpm+1)', fill = 'Ligand') +
  scale_fill_viridis_d(option = 'turbo', begin = .1, end = .9)
